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Chiral spin-orbital liquids with nodal lines 
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Strongly correlated materials with strong spin-orbit coupling hold promise for realizing topolog¬ 
ical phases with fractionalized excitations. Here we propose a chiral spin-orbital liquid as a stable 
phase of a realistic model for heavy-element double perovskites. This spin liquid state has Majorana 
fermion excitations with a gapless spectrum characterized by nodal lines along the edges of the Bril- 
louin zone. We show that the nodal lines are topological defects of a non-Abelian Berry connection 
and that the system exhibits dispersing surface states. We discuss some experimental signatures of 
this state and compare them with properties of the spin liquid candidate Ba 2 YMo 06 . 


PACS numbers: 71.70.Ej, 75.10.Jm, 75.10.Kt 

Quantum spin liquids (QSLs) are Mott insulators in 
which quantum fluctuations prevent long-range magnetic 
order down to zero temperature [1]. They have re¬ 
ceived both experimental and theoretical attention due 
to predictions of unusual phenomena such as spin-gapped 
phases with topological order or gapless phases without 
spontaneous breaking of continuous symmetries [2]. In 
recent years the evidence for QSLs in nature has started 
to look more auspicious thanks to the discovery of new 
compounds that realize the Heisenberg model on frus¬ 
trated lattices [3] . While frustration is a desirable ingre¬ 
dient, seminal work by Kitaev [4] has demonstrated that 
bond-dependent exchange interactions may provide an¬ 
other route towards QSL ground states. The key idea is 
that a spin-1/2 model on the (bipartite) honeycomb lat¬ 
tice with judiciously chosen anisotropic interactions can 
be rewritten in terms of free Majorana fermions hopping 
in the background of a static Z 2 gauge held. The result 
is a QSL with exotic fractional excitations. The same 
idea has been applied to construct other exactly solvable 
models, including cases of higher spins [5-9]. 

From a broader perspective, Kitaev’s model is an in¬ 
stance of a quantum compass model [10-12]. Although 
Kitaev-type exactly solvable models are artificial, the 
kind of anisotropic interactions they presuppose arises 
naturally in Mott insulators with orbital degeneracy and 
strong spin-orbit coupling [13, 14]. There is recent evi¬ 
dence that bond-dependent interactions are dominant in 
Na 2 lr 03 [15] . While this compound is in a zigzag-ordered 
phase at low temperatures, the prospect of finding QSLs 
in compass models suggests inspecting other families of 
heavy-element transition metal oxides [16-18]. 

All the conditions leading to quantum compass mod¬ 
els can be found in Mott-insulating rock-salt-ordered 
double perovskites [19]. Given the chemical formula 
A 2 BB' 06 , particularly interesting properties are found 
in compounds where the B' magnetic ions have a 4d^ 


or 5d^ configuration. These ions are arranged in a face- 
centered-cubic (fee) lattice, which, unlike the honeycomb 
lattice, is geometrically frustrated. The magnetic prop¬ 
erties within this family are diverse [20-23] , but the ma¬ 
terial that stands out is Ba 2 YMo 06 [24-27]. Despite a 
Gurie-Weiss temperature 0cw ~ —160 K [24] , several ex¬ 
periments point to the absence of long-range order down 
to T ^ 2 K. Moreover, there is no sign of structural tran¬ 
sitions, implying that the lattice remains cubic at low 
temperatures. Thus, the degeneracy of the t 2 g orbitals is 
only partially lifted by the spin-orbit coupling, leading to 
a low-lying j = 3/2 quadruplet [24] . The effective model 
contains bond-dependent interactions between nearest- 
neighbor j = 3/2 spins and is closely related to T-matrix 
generalizations of Kitaev’s model [5, 7[. Remarkably, the 
analysis in [19] revealed that, when antiferromagnetic 
exchange is dominant, ordered phases become unstable 
against quantum fluctuations, making this an interesting 
point to look for QSLs. 

In this Letter we investigate a QSL in a realistic model 
for double perovskites with strong spin-orbit coupling. 
Using a representation of j = 3/2 operators in terms of 
six Majorana fermions, we start by showing that a hidden 
SU (2) symmetry of the Hamiltonain becomes an explicit 
SO(3) symmetry for three of these fermions, whereas 
the other three exhibit a compass-model-type Z 3 sym¬ 
metry. As the model is not exactly solvable, we proceed 
with a mean-field approach and propose a spin liquid 
ansatz that preserves the SO(3) and Z 3 symmetries. The 
ansatz breaks inversion and time reversal symmetry, thus 
describing a chiral spin liquid. Most interestingly, we 
And that the excitation spectrum is gapless along nodal 
lines which are vortices of a Berry connection in momen¬ 
tum space. This feature makes this chiral spin liquid a 
strongly correlated analogue of line-node semimetals and 
superconductors discussed in the context of topological 
phases of weakly interacting electrons [28-31] and pho- 
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tonic crystals [32] . Going beyond the mean-field level, we 
use variational Monte Carlo (VMC) techniques [33-35] to 
show that our spin liquid state yields a remarkably low 
energy and should be regarded as a competitive candidate 
for the ground state of the spin-orbital model. Finally, we 
argue that the vanishing density of states at low energies 
predicted by our theory can account for some unusual 
properties observed in Ba 2 YMo 06 . 

The spin-orbital model for cubic double perovskites 
with d^ electronic configuration is given by [19] 

H — J ^ ^ f j ^ ^ ^ li ' Sj. 

a,{ij)£Q. i 

( 1 ) 

Here J > 0 is the antiferromagnetic exchange between 
nearest-neighbor spins and A > 0 is the atomic spin-orbit 
coupling. The index a labels both planes (XY, YZ or XZ) 
and t 2 g orbitals [d^y, dy^ or d^z) [36] . The operators rii^a 
and Si_Q, describe the number and the spin of electrons 
occupying the a orbital on site i, with the constraint 
X)a effective I = 1 orbital angular 

momentum of the t 2 g states [37] . The total spin on site 
i is Sj — 

In the regime A 3> J, spin and orbital operators can be 
projected into the low-energy subspace of total angular 
momentum j = 3/2 [19[. The projected Hamiltonian 
H = 7^3/2777^3/2, where V 3/2 is the projector, contains 
multipolar interactions in terms of = 1^ -l-S^. Our first 
step is to introduce operators s and r at each site as 

i(_r23,ri3^ri2), t = ( 2 ) 

The notation refers to five Dirac T matrices given explic¬ 
itly by 

pl = cr^ (g) ay, (g) T^ = a!' (g) 1, 

T"* = cr“ (g 1 , ( 3 ) 

where (t“, a S {x, y, z}, are Pauli matrices, and 10 matri¬ 
ces = [r'^,r‘^] /{2i) [5, 38] . The components of s and 
T satisfy the SU(2) algebra [s“,s^] = = 

ieab^rC, and [s“,T^j = 0. The relation between the ba¬ 
sis of and the basis |s^,t^), with s^,r^ G {+,—}, is 

|j^=±|) = |T,+),|J^=±i) = |±,-). 

In the new representation the projected Hamiltonian 
assumes a relatively simple form: 

^ ^ E - i) (1 - - 2^“)- (4) 

{i,j)ea 

where r“ are given by = t*, = 

7(—±a few comments are in order. First, 
Eq. (4) has the familiar form of a Kugel-Khomskii model 
[39, 40] . However, here the Kugel-Khomskii coupling in¬ 
volves pseudospins s and pseudo-orbitals t defined in the 
j = 3/2 subspace, where the original spins and orbitals 


are highly entangled. Second, the Hamiltonian commutes 
with Stot = This is a manifestation of the hidden 

global SU(2) symmetry discussed in [19] . This continuous 
symmetry is unexpected, given that spin-orbit coupling 
breaks the conservation of Jtot = J*) but appears in 

related models for t 2 g orbitals [41] and at special points 
of the Kitaev-Heisenberg model [42] . Finally, the pseudo¬ 
orbital coupling has the form of a 120 ° compass model 
[12[. There is a Z 3 symmetry generated by C /3 = 
followed by a C 3 rotation of the a planes. 

In analogy with the spin liquid in Kitaev’s model [4], 
we now introduce a Majorana parton representation for 
the generators of SU(4) (i.e. the basis of j = 3/2 opera¬ 
tors). We write s and t operators as [43-48] 


s“ = --e^bc„b c 

bj q ° h h ’ 


a ^ ^abcnbnc 


( 5 ) 


The components of rjj and 6 j are Majorana fermions that 
obey { 7 “, 7 /} = 26jiS°‘b, where 7 G {77, 9}. As the signs of 
the fermions can be changed (77 — —rj, 6 —>■ — 0 ) without 
affecting the physical operators, this representation bears 
a Z 2 redundancy. To eliminate the extra states, one can 
impose the local constraint [45] 

Dj = Dg = 1 Vj. ( 6 ) 


With this constraint we have s^rj* = —jVjd’j , and Hamil¬ 
tonian (4) becomes quartic in Majorana fermions: 


H 



HTihW, + vhWj + vhlvjWj + (* ^ j) 
+9:9fm ■ m , 


( 7 ) 


where N is the number of sites and 9°‘ and 0“ are defined 
by e^y = 01 , = i(_ 6 )i =p y/le^), and e^y = e^, 

= l(-03±v^0i). 

Hamiltonian (7) is invariant under global SO(3) rota¬ 
tions of the T 7 vector. The couplings involving the com¬ 
ponents of 0 have only a discrete symmetry, namely the 
octahedral point group symmetry Oh of the lattice. The 
latter contains the Z 3 that rotates 0“ and 0“ by 120° in 
the (0i,03) plane. In addition, H is invariant under time 
reversal T = where K denotes complex conju¬ 

gation. In terms of Majorana fermions, T = K 

Next, we perform a mean-field decoupling of Hamilto¬ 
nian (7). This is equivalent to neglecting fluctuations of 
the Z 2 gauge field and yields qualitatively correct results 
as long as the system is in a QSL phase with deconfined 
Majorana fermions [43] . Our choice of mean field ansatz 
is guided by the condition of preserving the SO(3) and Z 3 
symmetries. This restricts the set of nonzero parameters 
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Figure 1: (color online), (a) Two gauge-inequivalent hopping 
configurations with the same flux of the Z 2 gauge field on all 
faces of a tetrahedron, (b) Four-sublattice ansatz on the fee 
lattice. The sign of the outward flux alternates between edge¬ 
sharing tetrahedra (represented by different color fillings). 

ilil])- We obtain 

Hmy = + ^ H + wl)r^y rii 

+2,iuji9pf - iwfiOpf - 

+3u]i + - wfiVji] , (8) 

where iuji = {’q'-'nl), ivji = {OjOf), iwfi = ( 6 '“ 6 »“), and 
iw^i = {dfdf) play the role of imaginary hopping ampli¬ 
tudes. Note that the symmetry implies decoupling of 77 “ 
and 9^ fermions at the level of bilinear terms; yet, 9^ and 
9^ remain coupled. 

Seeking a translationally invariant ansatz, we set the 
order parameters to have uniform magnitude: Uij = uafj, 

C/gxy = with 

u,v,w,w to be determined by self-consistent equations, 
whereas the cr’s are chosen to be ±1 on each bond. Since 
e.g. Uij = —Uji, the choice of tjT is equivalent to a choice 
of bond orientation and determines the gauge-invariant 
flux through elementary plaquettes. Noticing that the 
fee lattice can be viewed as a network of edge-sharing 
tetrahedra, we obtain a symmetric ansatz by requiring 
that the Z 2 fluxes, e.g. x'jki — same 

on all faces of a given tetrahedron, with sites jkl on ev¬ 
ery triangle oriented counterclockwise with respect to an 
outward normal vector. This leads to the four-sublattice 
ansatz illustrated in Fig. 1. 

Let us discuss the symmetry of our ansatz. First, we 
note that the Z 2 gauge flux through triangles is odd un¬ 
der time reversal and is related to the spin chirality or¬ 
der parameter [49, 50]. The state also breaks inversion 
P; this can be seen from Fig. 1(b), since a mirror-plane 
reflection exchanges tetrahedra with opposite chiralities. 
Thus, our ansatz describes a chiral spin liquid with spon¬ 
taneous breaking of P and T. However, PT is still a 
symmetry. Similarly, a projective symmetry group anal¬ 
ysis [ 2 ] shows that broken rotational symmetries can be 
combined with the broken time reversal to restore an Oh 
point group symmetry, ensuring the orbital degeneracy 


Figure 2: (color online), (a) Brillouin zone of the cubic sub¬ 
lattice. (b) Bulk dispersion. The nodal lines along MR direc¬ 
tions cross at the quadratic band touching point R. 

assumed at the outset (see Supplemental Material [51]). 

Having fixed the ansatz, we can calculate the resulting 
spectrum of the Majorana fermions. For simplicity, first 
we focus on the mean-field Hamiltonian for 7 G { 77 “, 9^}, 
i.e., the flavors which are decoupled in Eq. ( 8 ). In this 
case the Hamiltonian can be written in the form 

^mf= E 7,^H(k)7k = |i| E 7,^(hk-S)7k, (9) 

kGjBZ kGjBZ 

where t = t(u,v,w,w) is the corresponding hopping 

amplitude in Eq. ( 8 ), 7 ^ = (tL- 7kB-7lc’7 Id) is 
a spinor with components labeled by sublattice index, 

hk = 4 ^cos ^ cos cos ^ cos cos ^ cos S = 
(-ri,-r3,ri3), and the sum is restricted to half Bril¬ 
louin zone since 7 _j. = 7 ^ [44] . As the components of S 
obey [S“, = 7 e“*'‘^E'^, the spectrum is given simply by 

e±(k) = ±|t||hk|. (10) 

The dispersion relation is illustrated in Fig. 2. There 
are two doubly degenerate bands [51]. Since {P, i^Mp} = 
0, the Hamiltonian has a chiral symmetry [52, 53] and the 
spectrum is symmetric between positive- and negative- 
energy states. The defining feature of the band struc¬ 
ture is the band touching along the edges of the Bril¬ 
louin zone. These are nodal lines parametrized, e.g., by 
k = ( 7 r, 7 r, fc^). Expanding k = (tt + p^, tt + py, k^), 
with Px,Py 1, we obtain the effective Hamiltonian 
on a plane perpendicular to a line node: P(k) « 
2|t| cos ^(pa;E^ + Py^^). The latter is formally equiv¬ 
alent to the Hamiltonian for graphene and yields lin¬ 
ear dispersion at low energies with fc^-dependent veloc¬ 
ity e±(k) « ±2|t| cos +Py- These nodal lines 

can be characterized as topological defects of an SU(2) 
Berry connection [54] in reciprocal space (see Supple¬ 
mental Material [51]). The three nodal lines related 
by C 3 symmetry cross at R = ( 7 r, 7 r, 7 r). Expanding 
k = {tt -\-pxiT^ PPz)-, we And that R is a quadratic 

band touching point [55, 56] with anisotropic dispersion 

£±(k) « \t\^plpl+plpl+plpl. 
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Figure 3: (color online). Surface-state spectrum projected in 
the Brillouin zone of the triangular lattice for (111) surface. 
The dispersion of the surface states corresponds to the thick 
red line between points M and K. 

Another feature of topologically nontrivial states of 
matter is the presence of protected surface states. We 
identify the surface states by calculating the spectrum 
for Hmf in n slab geometry with open boundary condi¬ 
tions in the (111) direction (Fig. 3). There appear two 
pairs of doubly degenerate bands separated from the con¬ 
tinuum, with dispersion terminating at the projections of 
the nodal lines. Remarkably, the positive-energy surface 
states are spatially separated from the negative-energy 
ones as their wave functions are localized at opposite sur¬ 
faces (which surface depends on the sign of the hopping 
parameter). This is a direct manifestation of the break¬ 
ing of inversion symmetry. 

We calculate the ground state energy Egg at mean-field 
level by solving the self-consistent equations that deter¬ 
mine the order parameters in Eq. (8). For this purpose 
we had to diagonalize the Hamiltonian for coupled 9^ and 
0^ fermions in Eq. (8) and found that the spectrum again 
displays nodal lines [51]. We obtain /NJ « —0.248. 
A better estimate of Egg can be obtained by implement¬ 
ing a Gutzwiller projection according to Eq. (6) using 
VMC [33, 45[. Considering a restrictive form of the wave 
function which neglects variations in the population of 
the fermionic flavors (see Supplemental Material [51[), 
we obtain E'^J^^/NJ = —0.40(1). This energy is al¬ 
ready comparable to that of the best variational state 
identified in Ref. [19[, namely a valence bond solid with 
Eybs/^J ~ —0.417. We expect the spin liquid to be 
stable since small fluctuations of the Z 2 gauge field only 
induce weak short-range interactions [2[ , which are irrel¬ 
evant in the renormalization group sense for topological 
semimetals with point or line band touching in three di¬ 
mensions [56[. 

The low-temperature thermodynamic properties are 
governed by the density of states p(e) oc of the Majo- 
rana fermions, which is due to the quadratic band touch¬ 
ing point. It follows that the QSL has heat capacity 
C oc magnetic susceptibility x oe and thermal 
conductivity k oc for ksT <C J. Another important 
property is the correlation function G(r) = (Jj • Jj+r). 


We And that G(r) vanishes when r connects sites on the 
same sublattice. For vectors connecting different sub¬ 
lattices along (100) directions in the form r = ^ -|- re, 
where 5 e {(i, §,0), (i,0, \), (0, i, |)} and e S {x,y,z}, 
the correlation decays at large distances as G(r) ~ l/r'^. 
This power-law decay coincides with the result for a Dirac 
point in two dimensions [45[. 

Finally, we address the comparison with available 
experimental results for the spin liquid candidate 
Ba2YMo06. Aharen et al. [25[ observed that both the 
heat capacity and the magnetic susceptibility vanishes at 
low temperatures and have attributed this behavior to a 
gapped collective spin singlet, de Vries et al. [24[ pro¬ 
posed a picture of a valence bond glass, but noted that 
the muon spin relaxation is comparable to that of QSLs 
[27[. Here we propose that an alternative explanation 
for the vanishing heat capacity and susceptibility at low 
temperatures is the vanishing density of states of our gap¬ 
less spin-orbital liquid with nodal lines. A comprehensive 
study of the properties of this QSL in comparison with 
experimental results will be presented elsewhere [57[. 

To summarize, we have studied a realistic model for 
double perovskites in the regime of strong spin-orbit 
coupling. We proposed a new spin liquid ansatz that 
gives rise to nodal lines in the spectrum of Majorana 
fermions. We argued that some experimental results for 
Ba 2 YMo 06 can be interpreted in terms of the vanish¬ 
ing density of states predicted by our theory. We hope 
this work will stimulate the search for strongly correlated 
materials hosting fractional excitations with nontrivial 
momentum-space topology [48, 58[. 
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Supplemental Material 
1. Symmetry of the spin liquid ansatz 

Time reversal T = = K ) acts on peudospins and pseudo-orbitals as 


T SjT — Sj, 


rp-l^X,Zrp _ 


T-^ryT = -ry. 


( 11 ) 


In the Majorana fermion representation for s and t this transformation can be implemented hy T = K n, 6»]6»3. This 
is equivalent to complex conjugation combined with the operation —>■ — and —> —dj- Thus, if we focus on the 
decoupled flavors 7 G { 77 “, 6 ^}, we can take time reversal to be represented simply by complex conjugation. 

Let jji denote a Majorana fermion on site j belonging to sublattice i = A,B,C,D. The operators 7 kf in momentum 
space are defined by 


Iji = 


= \ n 




TV 

keiBZ 




7kt = 


- VvE 


Ijie 


ik-R. 


( 12 ) 

(13) 




and 7 k^ are normalized such that { 7 k £5 7k'£'} “ For each flavor of Majorana fermion we combine the four 

sublattice modes into a single “spinor” 7 k = (7kA, 7kBj 7kC)Tko)*- 

In momentum space, time reversal takes k —)■ —k. Up to a hopping amplitude (determined by self-consistent equa¬ 
tions, see next section), the mean-field Hamiltonian for a decoupled flavor is of the form Hmf = X^ke^BZ Tk^(l^)7k 
with 


'H(k) = i 


( 0 f{K,ky) f{ky,k^) f{k^,k^) \ 

• 0 


-f(kx,ky) 0 -f{kx,kz) f{ky,kx) 


-f{ky,kz) fikx,k^) 


0 


-f{kx,ky) 


\-f{kx,kz) -f{ky,kz) f{kx,ky) 0 

where f{ka, kb) = 4cos(fca/2) cos(A:t,/2). Notice the factor of i. It follows that 

T-^H{k)T = •H*(-k) = -•H(k). 


J 


(14) 


(15) 


We define inversion P as the reflection by the mirror plane that exchanges A and C sublattices (plane perpendicular 
to Syz = (0, |)). In momentum space, P : kx kx,ky —>■ —kz,kz —>■ —ky. In addition, we have the action in the 

internal (sublattice) space given by the matrix (with determinant - 1 ) 


P = 


/O 0 1 0 \ 
0 10 0 
10 0 0 
\ 0 0 0 1 / 


(16) 


We also define the Z 2 gauge transformation that changes the sign of fermions on the B sublattice: 


Gp = 


/I 0 0 0 \ 
0-100 
0 0 10 
\ 0 0 0 1 / 


It is easy to check that inversion anticommutes with the mean-field Hamiltonian: 

{PGp)-^n{Pk)PGp = -n{k). 

It follows that the combined transformation PTGp is a symmetry of the Hamiltonian: 

iPTGp)-^n{Pk)PTGp = n{k). 


(17) 

(18) 

( 19 ) 
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The C 3 rotation about a (111) axis that leaves an A site invariant is represented by 


and the rotation in momentum space takes 
not required; we obtain immediately that 


M 0 0 0 \ 

0 0 0 1 
0 10 0 ’ 

V 0 0 1 0 / 

kz,ky —)■ kx,kz —> ky. In this case a gauge transformation is 


C3-lH(i?3k)C'3 = H(k). 

The C 2 rotation along the 2 axis that exchanges AoB, Co D is represented by 

/O 1 0 0 \ 

10 0 0 

*"2 “ 0 0 0 1 ’ 

\ 0 0 1 0 / 


and in momentum space i ?2 : k^ 


—h h 
•) 


—ky. We need to combine the C 2 rotation with the gauge transformation 


10 0 0 
0-100 
0 0 10 

0 0 0 -1 


We then have 


(C2G2)-iH(ii’2k)(C'2G2) =H(k). (24) 

Translation by 6 xy = ( 5 , ^,0), which we denote by T^y, has the same effect of exchanging sublattices as the above 
C 2 rotation. Thus, conjugation by TxyG 2 , together with Rj —>■ Rj + S^y in real space, is also a symmetry of the 
Hamiltonian (and likewise for the equivalent translations in yz and xz planes). 

Now consider a C4 rotation along the z axis going through an A site, which exchanges C and D sublattices. This 
can be represented in sublattice space by 

/ 1 0 0 0 \ 

G. = 0 10 0 r25l 

G4 0 0 0 -1 ■ 

\ 0 0 1 0 / 


In momentum space, i ?4 : k^ 


—kx- Combining with the gauge transformation: 


we find 


10 0 0 
0-100 
0 0-10 
0 0 0 1 


(G4G4)-^-H(R4k)(G4G4) = -H{k). (27) 

Thus, like P and T, the C4 rotation inverts the chirality of the ansatz. It is then easy to see that C4G4T is a symmetry 
of the Hamiltonian. 

The C4 rotation can be used to construct a symmetry transformation that accounts for the twofold degeneracy of 
the Majorana fermion bands. 

It is also interesting to consider the C4 rotation that exchanges A and B sublattices, given by 

/O -1 0 0 \ 

r' - 10 0 0 

G 4 - 0 0 10 • 

I 0 0 0 1 / 



If we define this to be a rotation around z axis in the opposite direction than the one in Eq. (25), the transformation 
in momentum space is R!^ = i.e., R'^ : kx ^ It is easy to check that the composition M = CiC'^ 

commutes with the Hamiltonian, M~^'H{'k)M = 'H(k), and obeys = —1. Thus, we can block diagonalize 'H(k) 
by sectors labeled by the eigenvalue ±z of the matrix M: 

= (29) 

where cr = (—cr*', cr®) and Um is the unitary matrix that diagonalizes M. It is then clear that the spectrum of 
■H(k) is twofold degenerate with eigenvalues ±|hk|. Two degenerate states can be distinguished by the eigenvalue ±1 
of the Hermitean matrix iM (which is analogous to the chirality of Weyl fermions in the massless Dirac equation). 

In summary, the chiral spin-orbital liquid ansatz lowers the symmetry of the Hamiltonian from OhxZ 2 (where Z 2 is 
time reversal) to Oh (where the new group contains combinations of broken point group symmetries with the broken 
time reversal). 


2. Berry connection 

The nodal lines can be characterized as topological defects of a Berry connection in reciprocal space. In our case, 
the Berry connection has to be non-Abelian due to the double degeneracy of the bands. Away from the nodal lines, 
we define the SU(2) connection 


^mn(k) = i(^/’m(k)|(9fe>„(k)), (30) 

where |'!/'m(k)) and |■^/;„(k)), m,n S {1,2}, are degenerate eigenstates of 'H(k) (say with energy e_|_(k)) chosen so as 
to obey ('!/'m(k)|'0„(k)) = dmn and to diagonalize A^(k). The generalized Berry phase is the Wilson loop 

U = V exp[—i dfcaA“(k)], (31) 

where V denotes path ordering. The calculation of U is simplified if we consider a path around the line node 
parametrized by k « (tt -|- e cos a,7r + e sin a, k^), with a € [—t^, t^)- For infinitesimal radius e <C 1, we obtain 

=-^sinaa2'-kC>(e°), ^ cos a a*'-k C>(e°), (32) 

which is precisely the singular e dependence of a vortex line. We then find C/ = — 1, equivalent to a tt Berry phase. 


3. Solving the mean-field Hamiltonian 

In this section we outline the steps required to diagonalize the mean-field Hamiltonian and calculate the ground 
state energy. 

Using the mode expansion Eq. (12), we can rewrite the various hopping terms for Majorana fermions in terms of 
operators in reciprocal space. For instance, 

* 7jA7/b = 2z Y ^k(7L7kB-7kB7kA)- (33) 

0,0 6XY keiBZ 

where h\^ = 4cos(fca;/2) cos{ky/2) is the first component of hk. The mean-field Hamiltonian becomes 

3 

{2u + w)Y ^i(k)7S - w {0iyni{k)6l + (0k)^ H 2 (k) 0 k 

a—1 

Y + 3w“u,j - w^yVij) , (34) 

{i,j)ecy. 

where 'Hi(k) is the 4x4 matrix in Eq. (14), 0k = (^^kAi ■ • ■ ^ ^ko)* I® eight-component spinor that combines 6 ^ 
and 9^ fermions and 7^2 (k) is an 8 x 8 matrix to be specified below. 


4 E 


18 
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First consider the fermions ( G {r]°‘,0^}, whose spectrum is determined by 'Hi(k). Let be the unitary transfor¬ 
mation that diagonalizes 'Hi(k): 

t/^Hi(k)C/k = Ai(k), (35) 

where Ai(k) = diag{—|hk|, — |hk|, |hk|, |hk|} is a diagonal matrix. The operators that annihilate fermions in eigen¬ 
states of 'Hi(k) are 

7kA = ^ {U^)\nki, (36) 

e=A,...,B 

4 

7k£ = (37) 

A=1 


where A = 1,..., 4 is the band index. The mean-field ground state |GS) is the state in which all single-fermion states 
with negative energy are occupied. This leads to the self-consistent equation for expectation values, e.g. 


(7j.A7j+5.„.B) 


^ ^ [e*'^'"“”(7lA7kB)+e (7kA7i^B)] 

keiBZ 


= ilm 


^ E E(^')>AC^BAe*‘^'*-(7L\A) 

keiBZ A 



V [ d^k{U^)xAUBX 


(38) 


where in the last line the sum is over bands with negative energy and we took the thermodynamic limit to replace 
Sk I (corresponding to A^/4 states in the Brillouin zone of the cubic sublattice). 

Since ?7i(k) determines the spectrum of 77 “ and 6 ^ fermions, the self-consistent equations for Uij = —i{r]iVj) ^-bcI 
Vij = —i{0i0'j) are the same up to an overall minus sign, depending on the relative sign of the hopping amplitudes 
2u + w and —w in Eq. (34) We then have the constraint |m| = |u|, but must analyze two possibilities, namely u = v and 
u = —V. Without loss of generality (by choosing one of the two degenerate ground states with opposite chiralities), 
we can set u > 0. Numerical evaluation of the integral in Eq. (38) then yields u « 0.258. 

The relation between u and v determines the 8 x 8 matrix for 0k. For v = u, we obtain 


H2(k) =uhk-S', 

where S' = {2K^ 0 , -2Ky 0 1 , -2K^ 0 a^), with 


/ 0 

2 

0 

73 \ 


/ 

0 

-2 

0 

73 \ 



\ 

2 

0 

0 

73 

73 

0 

0 

0 

_ 

, Ry = - 
’ 2 


2 

0 

0 

73 

-73 

0 

0 

0 

, = - 
2 

-1 

-3 

Vv^ 

0 

0 

0 


1 - 

-73 

0 

0 

0 


[ 

3/ 


(39) 


(40) 


The components of the matrix vector K satisfy the SU(2) algebra. We then obtain the spectrum of 7^2 (k) and use 
it to solve the self-consistent equations for Wij = —i{9f0^) and Wij = —i{9f9^) analogous to Eq. (38). In this case 
of u = w we find a self-consistent solution with w « —0.081 and w « 0.317. Having fixed the order parameters, we 
obtain the mean-field ground states energy £^¥{ 1 ^ = u) = (Hmf) ~ —0.244A^J. 

For V = —u we obtain 


= uhk • S", (41) 

where S" = ((2 — tr®*') 0 S^, (2 — 0 S^, (2 — cr^^) 0 E^), with S the matrix vector defined in the main text, and 

the 2x2 matrices cr“ given by ± a/Sct^). In this case we find a self-consistent solution 

with w « 0.161 and w « 0.318. The ground state energy is Eupiy = —u) ~ —0.248WJ, slightly lower than the result 
for u = u. This is the value quoted in the main text. We note that the small difference between the two energies may 
change beyond the mean-field level. However, we have verified that both solutions give rise to a spectrum with nodal 
lines along MR directions, qualitatively similar to the spectrum for 77 “ and 9^ fermions. Therefore, the properties 
derived from the low-energy density of states / 3 (e) oc y/e are generic. 
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4. Variational Monte Carlo 

To check the viability of the proposed chiral spin-orbital liquid beyond the mean-field level, we now enforce the 
local constraint exactly by considering a Gutzwiller projection of the mean-field wave function [33] by means of a 
variational Monte Carlo calculation [34] . 

We begin by rewriting the Majorana fermions in terms of three Dirac fermions, closely following the representation 
used in Ref. [45] 


4 = ^ + iOl) . 4 = \ (4 + i4) ^ fl = \ (4 + ■ 


(42) 


In terms of this representation, the local constraint is now translated into the fact that a given site may either have 
no Dirac fermions, a state denoted by |0), or two Dirac fermions, in which case there are three possible states at each 
site defined as 


\xi) = c\4 |0), \yi) = 4fl |0), \zi) = f\c\ |0). 


(43) 


Given a real space configuration specified by the locations of the doubly occupied sites, X = {xi} , Y = {?/j} , Z = {zi}, 
the wave function assigns an amplitude {{xi} ,{yj} ,{zm}) to it. Notice that the locations of the |0) states are 
automatically specified. We point out that the local constraint only fixes the parity of the number of fermions but 
not the number itself. Moreover, our Hamiltonian contains terms which not only create/annihilate two particles, e.g. 
\xi) ^ jo), but also terms which preserve the total number of fermions while changing the number of each of the three 
individual fermionic flavors, e.g. \xi) ^ \yi)- Although it is possible to write down a projected wave functions with 
varying particle number [35], we refrain from doing so in this work for the sake of computational simplicity. Instead, 
we consider a restrictive form for the ground state wave function: each one of the four states is equally distributed 
over iV/4 distinct sites, with {xi} = {^xi,X 2 , ■ ■ ■ ,xn/ 4 ^, etc. 

The mean-field Hamiltonian in the main text may be rewritten, in terms of the three Dirac fermions in Eq. (42), as 


T~Lmf = — 


NJ J 
IT 3(5 


+ SwjiUji - wfiVji) + i (Auji + 2wfi) 4di 




+ i {2uji + w'^i - wfi) c]ci + i {2uji + wfi + wfi) cjc) -k i {3uji - Vji) f] fi + h.c. 


36 

J 


Ui)eXY 


J 


X * 


{ji)eYZ 




)/]//+h.c.' 


X * 


maxz 


„-2i7r/3\ ft 


- v,ie I 


f +h.c. 


(44) 


The three fermion flavors are decoupled and we may thus write the mean-field wave function as a product of three 
Slater determinants. Thus, after the Gutzwiller projection, we obtain 


Avj} ,Um}) = ,{yj}) ■ {{Vj} ,{zm}) ■ ^c{{zm} ,{xi}). (45) 

The d-fermion sector of the Hamiltonian in Eq. (44) corresponds to free fermions and thus their mean-field ground 
state is obtained by filling up the states with negative energy. For the c and /-fermions we have a BGS-like Hamiltonian 
instead and their ground state is given by the vacuum of their respective Bogoliubov quasiparticles [33] . The different 
status of the d fermion is expected from symmetry: the hidden global SU(2) symmetry of the original Hamiltonian 
implies the global U(l) symmetry corresponding to the conservation of the total number of Dirac fermions defined 
by a combination of y Majorana fermions. On the other hand, there is no continuous symmetry associated with 
0 Majorana fermions; as a result, the total number of c and / fermions is not conserved. The need to work with 
BGS-type wave functions in our case should be contrasted with the case of SU(4) symmetric models [45], where the 
SU(4) symmetry implies the conservation of the numbers of all three flavors of Dirac fermions. 

After constructing the mean-field wave function, we then implemented a variational Monte Garlo calculation of the 
Gutzwiller-projected ground state energy . We started by generating an initial state in which we populate N/4 

randomly chosen sites with the x-state ({x^}), then 7V/4 of the remaining sites with the y-state ({j/i}), and finally 
A^/4 of the further remaining sites with the z-state {{zi}). Our Monte Garlo moves consist in exchanging random 
pairs of sites containing distinct states. We allow for moves involving widely separated sites — and which would not 
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be connected by the Hamiltonian — because this improves the sampling over the space of configurations. We accept 
or reject these moves according to the usual Metropolis algorithm. After N such exchange attempts, we are said 
to have performed one Monte Carlo sweep and after every sweep we compute . fVwarm sweeps are performed 

before measurements of physical quantities for “thermalization”. Averages are then performed over A^sweep sweeps. 
We typically considered Awarm = Asweep 10^. The results were obtained for lattices of size A = with L = 4, 
6 , and 8. We find that the change in the ground state energy with A is smaller than the Monte Carlo error bars for 
the system sizes considered here. Thus, we quote the results for L = 8 as the converged ones. 

We computed the ground state energy for the two sets of mean-field parameters quoted in this supplemental material. 
For u = V = 0.258, w = 0.317, and w = —0.081 we obtain = —0.39 (1) AJ. As for u = —v = 0.258, w = 0.318, 

and w = 0.161 we obtain = —0.40 (1) AJ. Clearly, the Gutzwiller projection decreases significantly the mean- 

field energy down to values which are already comparable to that obtained, for instance, for a valence-bond covering 
of the lattice (Avbs = —0.417AJ) [19], thus showing that the proposed chiral spin-orbital liquid is a competitive 
ground state candidate. 

In light of this favorable energy of our proposed ansatz, we conclude by pointing out two important restrictions in 

our variational Monte Carlo calculation that, once lifted, should further decrease the value of the ground state energy 
jpVMC. 

1. For the quoted values of , we considered the optimal values for the mean-field amplitudes v, w and w 

obtained before the Gutzwiller projection, i.e., at the mean-field level; 

2. The restrictive form of the considered wave function neglects variations both in the populations of the fermionic 
flavors and in the total number of fermions. 

We stress that these restrictions were important for this first calculation beyond mean field due to the complexity of 
the chiral spin-orbital liquid ansatz considered here. We leave a more detailed investigation, together with a more 
precise estimate for the variational energy of our spin liquid ansatz, for future work. 


